#ifndef PCPS_MC_ENERGY_HEADER
#define PCPS_MC_ENERGY_HEADER

#include <src/pcps.h>
#include <src/input.h>
#include <src/hamiltonian.h>
#include <src/wavefunction.h>
#include <src/accumulator_energy.h>
#include <src/stochastic_optimization.h>

namespace pcps {

  //--------------------------------------------------------------------------------
  // pcps::mc_energy -- Computes and prints the energy of the cps wavefunction
  //                    via MC sampling.
  //--------------------------------------------------------------------------------

  template <class S, class PART, class REF>
  inline S mc_energy(const pcps::Input & userinp, const pcps::Hamiltonian<S> & ham, pcps::Wavefunction<S> & wfn) {

    // divide the sample length among the different threads
    const long int sample_length = 1 + userinp.sample_length() / ( MPI::COMM_WORLD.Get_size() * userinp.nthreads() );

    // define a short name for a shared pointer to a process accumulator
    typedef boost::shared_ptr<pcps::ProcessAccum<S, PART, REF> > ProcAccPtr;

    // initialize a vector of accumulators and place the energy accumulator in it
    std::vector<ProcAccPtr> accumulators;
    { ProcAccPtr temp( new pcps::EnergyProcessAccum<S, PART, REF>(userinp, ham, wfn, sample_length) ); accumulators.push_back(temp); }

    // compute the energy via stochastic sampling
    pcps::stochastic_iterations<S, PART, REF>(userinp, ham, wfn, accumulators, 1, 2, 2, sample_length, userinp.burn_in_length());

    // return the energy
    S retval = ( dynamic_cast<pcps::EnergyProcessAccum<S, PART, REF> *>(accumulators.at(0).get()) )->total_energy();
    return retval;

  }

} // end namespace pcps

#endif
